Optimized Multichannel Quantum Defect Theory for cold molecular collisions 
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Multichannel quantum defect theory (MQDT) can provide an efficient alternative to full coupled- 
channel calculations for low-energy molecular collisions. However, the efficiency relies on interpola- 
tion of the Y matrix that encapsulates the short-range dynamics, and there are poles in Y that may 
prevent interpolation over the range of energies of interest for cold molecular collisions. We show 
how the phases of the MQDT reference functions may be chosen so as to remove such poles from the 
vicinity of a reference energy and dramatically increase the range of interpolation. For the test case 
of Mg+NH, the resulting optimized Y matrix may be interpolated smoothly over an energy range 
of several Kelvin and a magnetic field range of over 1000 G. Calculations at additional energies and 
fields can then be performed at a computational cost that is proportional to the number of channels 
N and not to JV 3 . 



I. INTRODUCTION 

Samples of cold and ultracold molecules have unique 
properties that are likely to have applications in many 
diverse areas. These include high-precision measurement 
[H H] j quantum information processing [3] and quantum 
simulation [4]. There is also great interest in the devel- 
opment of controlled ultracold chemistry [5]. 

Atomic and molecular interactions and collisions are 
crucial to the production and properties of cold and ultra- 
cold molecules. However, quantum-mechanical molecular 
collision calculations can be computationally extremely 
expensive. Such calculations are usually carried out using 
the coupled-channel method, in which the wavefunction 
is expanded 

JV 

V(r,T)=r- 1 J2 i Pi(T)iii(r). (1) 

i=l 

Here the N functions form a basis set for the mo- 

tion in all coordinates, r, except the intermolecular dis- 
tance, r, and ipi(r) is the radial wavefunction in channel 
i. Substituting this expansion into the time- independent 
Schrodinger equation and projecting onto the basis func- 
tion ifj (r) yields a set of N coupled differential equations. 
The properties of completed collisions are describe by the 
scattering matrix S, which is obtained by matching the 
functions ipi(r) to free-particle wavefunctions (Ricatti- 
Bessel functions) at long range [6j. In the full coupled- 
channel method, explicit solution of the coupled equa- 
tions takes a time proportional to A 3 . 

The problems encountered in cold molecular collisions 
often require very large number of channels. Atom- 
molecule and molecule-molecule interaction potentials 
can be strongly anisotropic, requiring large basis sets of 
rotational functions for convergence. In addition, calcu- 
lations are often required in an applied field, where the 
total angular momentum J is no longer a good quantum 



number. Because of this, the large sets of coupled equa- 
tions cannot be factorized into smaller blocks for each J 
as is possible in field-free scattering [7]. Furthermore, at 
the very low collision energies of interest, small splittings 
between molecular energy levels have important conse- 
quences. Effects such as tunneling [8j and nuclear hyper- 
fine splitting pillOj each multiply the number of channels. 

In cold collision studies, the scattering S matrix is of- 
ten a fast function of collision energy E and magnetic 
field B, with extensive structure due to scattering res- 
onances and discontinuous behavior at threshold. Cal- 
culations are often required over a fine grid of energies 
and/or applied electric and magnetic fields, and this fur- 
ther multiplies the computational expense. 

We have recently shown [IT] that Multichannel Quan- 
tum Defect Theory (MQDT) [PW7] provides an attrac- 
tive alternative to full coupled-channel calculations in 
these circumstances. MQDT attempts to represent the 
scattering properties in terms of a matrix Y{E,B) [T4l - 
ITT] that is a smooth function of E and B. If this can 
be achieved, the matrix can be obtained once and then 
used for calculations over a wide range of energies and 
fields, or obtained by interpolation from a few points. 
Once the matrix Y(E, B) has been obtained, the time 
required for calculations at additional energies and fields 
is only proportional to A, not A 3 . 

One problem with MQDT is that the Y matrix may 
have poles as a function of E and B, and these limit the 
range over which it can be interpolated. In cold molecular 
collision studies, calculations are typically needed over an 
energy range of order 1 K above threshold, and for mag- 
netic fields up to a few thousand gauss [18] . This con- 
trasts with the situation for collisions of ultracold atoms, 
where the energy range of interest is commonly a few /jK. 
and the fields are typically a few hundred gauss. 

In the present paper, we show how MQDT Y matrices 
can be defined to allow smooth interpolation over sub- 
stantial ranges of collision energy and applied field. This 
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will allow the use of MQDT to provide substantial sav- 
ings in computer time. 

II. THEORY 

A full description of MQDT has been given previously 
P2HI1]- We give here only a brief description, following 
ref. [TT], which is sufficient to describe the notation we 
use. 

MQDT defines the matrix Y(E, B) at a matching dis- 
tance r matc h at relatively short range. The iV-channel 
scattering problem at energy E is partitioned into iV 
open channels (with Ef° < E, where E°° is the threshold 
of channel i), N c weakly closed channels, and N s strongly 
closed channels. Strongly closed channels are those that 
make no significant contribution to the scattering dynam- 
ics at r > r matc h- 

The scattering dynamics beyond r matc h is accounted 
for using single-channel (uncoupled) calculations in a ba- 
sis set that diagonalizes the Hamiltonian at r = oo. 
The solution of the multichannel Schrodinger equation 
at r > r matc h is written in the matrix form 



* = f{r)+g(r)Y, 



(2) 



where / and g are diagonal matrices containing the func- 
tions fi and <7j, which are linearly independent solutions 
of a reference Schrodinger equation in each asymptotic 
channel i, 



+ uT(r) 



Mr) = 0, 



(3) 



and similarly for <?,(r). The reference potentials U\ (r) 
approach the true potential at long range, and /i is 
the reduced mass. They include the centrifugal terms 
h 2 Li{Li + l)/2/j,r 2 , where Li is the partial-wave quan- 
tum number for channel i. Y is an iV re f x iV re f matrix, 
where N rc{ = N + N C . 

In our approach Y is obtained numerically by 

matching the solutions of the coupled-channel equations 
to fi(r) and gt(r) at r ma t c h- The S matrix is then 
obtained from Y using Eqs. (21) to (23) of ref. [TT] . 
which require 3 QDT parameters Cj, tanA^ and £j in 
each open channel and a single QDT parameter tan Vi 



in each weakly closed channel. In the open channels the 
reference functions are asymptotically related to Ricatti- 
Bessel functions J.L 4 (r) and N L .(r) [5], 



-C'i tan A; 





cr 



cos^ sin^ 
-sin& cos& 



N 



Li 



(4) 

Here & is the asymptotic phase shift of the func- 
tion fi with respect to the Ricatti-Bessel function J^. 
The QDT parameter Q relates the amplitudes of the 
energy-normalized functions at long range to functions 
with Wentzel-Kramers-Brillouin (WKB) normalization 
at short range, while tan Xi describes the modification 
of the WKB phase due to threshold effects. Far from 
threshold, C, « 1 and tanA^ ss 0. In the weakly closed 
channels the reference functions are asymptotically 



cos Vi sin Vi 
- sin Vi cos Vi 



(5) 



where (f>i is the solution of ([3]) that decays exponentially 
at large r and 7$ is its linearly independent partner, which 
is exponentially growing. 

The absolute phases chosen for the reference functions 
fi and gi are arbitrary, and different choices produce dif- 
ferent Y matrices and different MQDT parameters. In 
particular, Eq. ^ shows that a pole in Y occurs when- 
ever the propagated multichannel wavefunction in any 
channel i has no contribution from the reference func- 
tion fi. However, all phase choices produce the same 
physical S matrix. We are therefore free to choose the 
phase in order to produce a Y matrix with advantageous 
characteristics. Here we show how the phase may be cho- 
sen to produce a Y matrix that is pole-free over a wide 
range of energy or magnetic field and can be interpolated 
smoothly. 

Rotating the reference functions fi and gi by an an- 
gle 9i gives a new set of linearly independent reference 
functions fi and <ji, 



cos tfj — sin ( 
sin 9i cos ( 



(6) 



These rotated reference functions define a new Y matrix 
and a new set of QDT parameters (C, tan A, £ and tan v) . 
Combining equations Q, ^ and ^ gives 



arctan 



sin £j (cos 8i + tan A^ sin 9i ) — cos sin ( 



c 



2 cos £i (cos 9i + tan A, sin + sin & sin 9i 



tanAi = — 



2Cf tan Xi cos 29, + [l + C? (tan 2 A, - l)] sin 29, 
2 (Cf cos 2 9i + sin 9{ [sin 9, + Cf tan Xi (2 cos 9i + tan Xi sin 9i)])' 



C, 



sin sin 9i 



. \ (cos&sinflj - C?sin&(cos0j + tan A, sin^)) 2 

■ (7j cos£j(cosPj + tanAjSinc'j) 1 . 1 H 2 

/ \ (sin£i sin^ + C 2 cos£i(cos&i + tan Ai sin^i)) 



(7) 
(8) 

(9) 
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(10) 



Far from threshold (£»1K), Eqs. Q to ([ToJ simplify 
to & = & — ^i) tan. A* ss 0, C 4 « 1 and P l = i/ ( - 0,. 
However, in the threshold region that is of interest in 
cold molecule studies, Eqs. ([7| to ( 10 1 must be evaluated 
explicitly. 



A. Basis sets and quantum numbers 

As a test case, we consider cold collisions between NH 
( 3 £~) and Mg atoms jTHj. This is the same system as 
considered in ref . [11] , but the present work uses a larger 
basis set which introduces more scattering resonances 
and denser poles in the Y matrix. 

The energy levels of NH in a magnetic field are most 
conveniently described using Hund's case (b), in which 
the molecular rotation n couples to the spin s to produce 
a total monomer angular momentum j. In zero field, 
each rotational level n is split into sublevels labeled by 
j. In a magnetic field, each sublevel splits further into 
2j + 1 levels labeled by rrij, the projection of j onto the 
axis defined by the field. For the n = levels that are 
of most interest for cold molecule studies, there is only 
a single zero-field level with j = 1 that splits into three 
components with rrij — +1, and — 1. 

The coupled equations are constructed in a partly cou- 
pled basis set \nsjmj)\LML,) , where L is the end-over-end 
rotational angular momentum of the Mg atom and the 
NH molecule about one another and Ml is its projec- 
tion on the axis defined by the magnetic field. Hyperfine 
structure is neglected. The matrix elements of the to- 
tal Hamiltonian in this basis set are given in ref. |20j . 
The only good quantum numbers during the collision are 
the parity p = (— 1)"+ L + 1 and the total projection quan- 
tum number M — rrij + Ml. The calculations in the 
present work are performed for p = — 1 and M = 1. This 
choice includes s-wave scattering of NH molecules in ini- 
tial state rrij — +1, which is magnetically trappable, to 
rrij = and —1, which are not. The present work uses a 
converged basis set including all functions up to n max = 6 
and L max = 8, as in ref. [19] . 

We label elements of Y and S by subscripts 
a,L,M L — > a',L',M' L , where a represents an eigen- 
state of free NH that may be approximately labeled by 
(n, s, j,rrij). However, the collisions considered in the 
present paper are all among the n = 0,j = 1 levels and 
so a is simply abbreviated to rrij. For diagonal elements 
we suppress the second set of labels. 



B. Numerical methods 

The coupled-channel calculations required for both 
MQDT and the full coupled-channel approach were car- 
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ried out using the MOLSCAT package [3T], as modified 
to handle collisions in magnetic fields [2U]. The cou- 
pled equations were solved numerically using the hybrid 
log-derivative propagator of Alexander and Manolopou- 
los [55], which uses a fixed-step-size log-derivative prop- 
agator in the short-range region (r mm < r < r m ;d) and a 
variable-step-size Airy propagator in the long-range re- 
gion (r mic j < r < r max ). The full coupled-channel calcu- 
lations used r m i n = 2.5 A, r m id = 50 A and r max = 250 A 
(where 1 A = 10 -10 m). MQDT requires coupled-channel 
calculations only from r m ; n to r ma tch (which is less than 
''mid), so only the fixed-step-size propagator was used in 
this case. 

The MQDT reference functions and quantum defect 
parameters were obtained as described in ref. [TT], us- 
ing the renormalized Numerov method [23] to solve the 
1-dimension Schrodinger equations for the reference po- 
tentials. The MQDT Y matrix was then obtained by 
matching to the log-derivative matrix extracted from the 
coupled-channel propagation at a distance r matc h. In this 
paper all MQDT calculations use the reference potential 



Ufir) = V (r) + 



ft 2 £,(L. t + l) 
2iir 2 



E° 



(11) 



where Vb(r) is the isotropic part of the interaction poten- 
tial. This reference potential has been shown to produce 
quantitatively accurate results when Y is reevaluated at 
each collision energy and magnetic field llj. However, 
such rccvaluation relinquishes most of the computational 
savings that MQDT is intended to achieve. 

The reference potential contains a hard wall at r = 
r, wa11 , so that U- cf (r) — oo for r < rf a11 . In the present 
paper we take rf a11 = 4.0 A. Figure [j] shows the refer- 
ence potentials for the lowest three rotational states. All 
channels with n > 2 were treated as strongly closed and 
thus not included in the MQDT part of the calculation, 
but were included in the log-derivative propagation. 



III. RESULTS AND DISCUSSION 

The top panel of Fig.[2]shows a single diagonal element 
of the Y matrix, Y"_ 1,8+3, as a function of the matching 
distance and energy, obtained with unrotated reference 
functions. YLx,8,+3 is a representative element of Y with 
poles at the same locations as the other elements, chosen 
to give a good visual representation of the pole structure. 
There are many poles visible, which prevent polynomial 
interpolation over energies of more that 0.5 K for any 
value of r matc h (and much less than this for some choices 
of r ma t c h). The energies of the poles become independent 
of r match at long range. 

The presence of low-energy poles in Y for some values 
of r match is a serious problem. For MQDT to be effi- 
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1(A) 



FIG. 1. The Vo reference potentials for Mg + NH. The first 
and second rotational excited state are also shown (n — 1,2). 
The hard wall at r = 4.0 A is shown as a vertical dashed line. 
The dot-dashed horizontal line corresponds to zero energy. 



cient, r ma t c h must be chosen without solving the coupled 
equations at many different energies. The calculations 
needed to produce contour plots such as those in Fig. [2] 
are feasible for a test case such as Mg+NH, but would 
be prohibitively expensive for a very large system. 

The center panel of Fig. [2] shows the same element of 
the Y matrix as a function of the matching distance and 
energy for reference functions rotated by 6i — tt/2. The 
poles are in quite different places, but once again there 
are many of them. The combination of the top and center 
panels demonstrates that, for any arbitrary choice of ro- 
tation angle, poles will appear in the Y matrix, prevent- 
ing simple interpolation for most choices of r ma t c h- This 
will be true in any MQDT problem with a large density 
of resonances. The contour plots do however show that 
the position of poles is strongly dependent on the rota- 
tion angle, even at large values of r matc h. This suggests 
that it will be possible to optimize the rotation angle in 
order to move the poles away from the energy range of 
interest. It is emphasized that the S matrices obtained 
from the Y matrices shown in the different panels of Fig. 
[2] are identical. 

We now consider how to rotate the reference functions 
to maximize the pole-free range over which Y can be 
interpolated. Ya as a function of 9i is given by 



Yu = tan(0* +Si), 



(12) 



where Si is the phase shift between the unrotated refer- 
ence function fa and the propagated multichannel wave- 
function in channel i. There is a pole in Yu when 
9i + Si = 7r and a zero when 9i + Si = 0. We thus set 
9° pt = —Si at one choice of r match , E and B, so that the 
propagated multichannel wavefunction and the reference 
wavefunctions are almost in phase and the resulting Y 
matrix in that region is pole-free. 

Because the channels are coupled, rotating the refer- 
ence functions in one channel affects the other elements 
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FIG. 2. (Color online) Contour plot of arctan y^/V for a rep- 
resentative diagonal Y matrix element, Y-1,8,+2, as a func- 
tion of energy and r matc h at B = 10 G. Top panel: obtained 
with unrotated reference functions (6i = 0). Center panel: 
obtained with reference functions rotated by 9% = 7r/2. Bot- 
tom panel: obtained with optimized reference functions with 
9i = 9° pt in all channels. The arctangent is show for clarity 
of plotting: it maps the real numbers, R, to the domain —n/2 
to 7r/2, thus allowing all magnitudes of Y matrix elements to 
be seen on a single plot. 
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of the Y matrix. In this work we loop over the channels 
sequentially, setting each diagonal element to in turn. 
By repeatedly looping over all channels, all the diagonal 
Y matrix elements are set to 0. For Mg+NH it was suffi- 
cient to loop over the channels twice. In a more strongly 
coupled system it is expected that this would need to 
be repeated more times. This approach allows a set of 
optimized 9i to be obtained from a single multichannel 
propagation. 

Rotated reference functions have previously been used 
to transform Y matrices in the study of atomic spec- 
tra [24-28 and atomic collisions [25]. Adjusting 6>; at 
each energy such that Ya = was shown to produce 
a weak energy dependence of off-diagonal Y matrix el- 
ements across thresholds [25J. However, this approach 
required propagating the full multichannel wavefunction 
many times at different energies, which is precisely what 
the present work tries to avoid. 

The bottom panel of Fig. [2j shows how the representa- 
tive clement Y°^ s +2 varies as a function of the match- 
ing distance and energy. All the 6i values are optimized 
as described above at E — 0.5 K and B = 10 G for 
each value of r matc h, but are not reoptimised at each 
energy. Comparison of this with the upper two panels 
shows the effectiveness of optimizing the reference func- 
tions. Without optimization, there were no choices of 
''match for which Y was pole-free and thus suitable for 
interpolation over the energy range of interest. After op- 
timization, Y opt is pole-free over a substantial range, of 
about 1 K, for any choice of r matc h < 8 A. For values 
of ^match < 6.5 A, Y° pt is pole-free over many Kelvin. 
Beyond 6.5 A, poles start to enter Y opt in the energy 
range of interest. Once the poles have settled at their 
asymptotic values at r m atch > 7.5 A, we find that pos- 
itive energies up to about 2 K are pole-free. However, 
at larger values of r ma tch the linearity of Y° pt over the 
pole- free region decreases. This is due to negative en- 
ergy poles in the Y matrix which our procedure cannot 
move significantly. There is one particularly bad choice 
of '"match at ~ 6.8 A, but provided this unlucky choice 
of '"match is avoided, Y opt can be interpolated smoothly 
over the positive energy range from to > 2 K for any 
choice of r matc h- 

Figure |3] compares diagonal T-matrix elements \Ta\ 2 
(where = Sij — Sij) obtained from full coupled-channel 
calculations with those from the MQDT method, with a 
matching distance of r matc h = 6.5 A, using reference func- 
tions optimized at 0.5 K. MQDT results were obtained 
both by recalculating the Y matrix at every energy and 
by interpolating Y opt linearly between two points sepa- 
rated by 1 K. The MQDT results with Y recalculated at 
each energy can scarcely be distinguished from the full 
coupled-channel results. The MQDT results obtained by 
interpolation are also very similar to the full coupled- 
channel results except around the resonance feature at 
E » 0.1 K. The interpolated result could of course be 
improved simply by performing coupled-channel calcula- 
tions to obtain Y opt at one or two extra energies across 




FIG. 3. (Color online) The squares of diagonal T-matrix el- 
ements T mj ,L,M L in the incoming channels for rrij = +1 and 
L — 0, 2 and 4 at B — 10 G, obtained from full coupled- 
channel calculations (solid, black) and MQDT with optimized 
reference functions for r matc h = 6.5 A, both with (dot-dash, 
blue) and without (dashed, red) interpolation. 




Bin 



FIG. 4. (Colour online) The asymptotic phase shift £j as 
a function of the rotation angle 9 for the incoming d-wave 
channel (+1,2,0). 



the range, to allow for a higher-order interpolation, or by 
using a linear interpolation over a smaller energy range. 

In this work we use 0i to rotate our short-range ref- 
erence functions /j and gi. In principle, we could rotate 
the reference functions by varying the asymptotic phase 
shifts £j instead of the short-range phases 0j. However 
Figure [4] shows why this is not desirable. Due to the 
highly nonlinear relationship between and 6i, obtain- 
ing the optimum rotation angle of the short-range refer- 
ence functions fi and by varying the angle £j would be 
laborious at very low collision energies. 
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FIG. 5. (Color online) Representative y opt matrix elements 
as a function of field at E = I mK. 



A. Magnetically tunable Feshbach resonances 

The effects of magnetic fields on cold molecular col- 
lisions are very important, since collisions can be con- 
trolled by taking advantage of magnetically tunable low- 
energy Feshbach resonances. We are therefore interested 
in how S matrix elements behave as a function of mag- 
netic field across Feshbach resonances. It is thus impor- 
tant that the Y matrix is weakly dependent on magnetic 
field in such regions. 

Figure [5] shows the diagonal elements of the optimized 
Y matrix as a function of magnetic field for Mg + NH 
collisions over the range from 10 G to 5000 G for a col- 
lision energy of 1 mK. This range of fields tunes across 
6 Feshbach resonances. The reference functions were op- 
timized at 10 G and 1 mK. The elements of Y opt are 
smoothly curved over the entire 5000 G range and could 
be well represented by a low-order polynomial. 

Figure [6] shows the comparison between optimized 
MQDT and full coupled-channel calculations for a se- 
lection of diagonal and off-diagonal T-matrix elements 
as the magnetic field is tuned at 1 mK. The reference 
functions were optimized at 10 G and 1 mK and MQDT 
results were obtained by linear interpolation of Y opt be- 
tween two points separated by 1000 G and by 5000 G. 
Interpolation over 1000 G gives resonance features that 
are in very good agreement with the full coupled-channel 
calculation to better than 1 G. Interpolation over 5000 G 
gives resonance features of the correct shape, with posi- 
tions that are still within about 10 G of the full coupled- 
channel results. The difference between the interpolated 
result and the full coupled-channel calculation is a result 
of both the choice of r ma t c h and the interpolation. The 
quality of the interpolation could be improved by con- 
sidering a few more fields across the range to allow for 
higher-order polynomial interpolation or by using linear 
interpolation over a smaller field range. 

Full MQDT calculations recalculating the Y matrix at 
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FIG. 6. (Color online) Squares of T-matrix elements at 1 mK 
as a function of field in the vicinity of a Feshbach resonance. 
Upper panel: diagonal elements; lower panel: off-diagonal 
elements. 



every magnetic field give resonance positions accurate to 
0.4 G. The remaining errors between the full coupled- 
channel calculations and the MQDT results will reduce 
with a larger choice of r matc h- As seen in the bottom 
panel of Fig. [2j the optimized Y matrices obtained at 
larger values of r matc h are still amenable to interpolation, 
though over a more restricted energy range. 



IV. CONCLUSIONS 

We have shown that Multichannel Quantum Defect 
Theory can provide an efficient computational method 
for low-energy molecular collisions as a function of both 
energy and magnetic field. In particular, we have shown 
how a disposable parameter of MQDT, the phase of the 
short-range reference functions, may be chosen to make 
the MQDT Y matrix smooth and pole-free over a wide 
range of energy and field. This smooth variation allows 
the Y matrix to be evaluated from coupled-channel cal- 
culations at a few values of the energy and field and then 
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to be obtained by interpolation at intermediate values. It 
is not necessary to repeat the expensive coupled-channel 
part of the calculation on a fine grid. 

The procedure developed here is to choose the phase of 
the reference functions in each channel so that the diago- 
nal Y matrix in each channel is zero at a reference energy 
and field. This ensures that there are no poles in the Y 
matrix, which would prevent smooth interpolation, close 
to the reference energy. Optimizing the phase in this 
way is very inexpensive, and once it is done the cost of 
calculations at additional energies and fields varies only 
linearly with the number of channels N, not as N 3 as 
for full coupled-channel calculations. MQDT with opti- 
mized Y matrices is thus a very promising alternative to 
full coupled-channel calculations for cold molecular col- 
lisions, particularly when fine scans over collision energy 
and magnetic field are required. 

The Y matrix is defined to encapsulate all the collision 
dynamics that occurs inside a matching distance r matc h, 
and the choice of this distance is important. There is a 
trade-off between the accuracy of the method and the size 
of the pole-free region of the optimized Y matrix. For 
large values of r matc h, resonant features may appear in 



the Y matrix and prevent simple interpolation over large 
ranges of energy and field. For smaller values of r matc h, 
optimizing the reference functions allows interpolation 
over many Kelvin, but the accuracy of MQDT is reduced 
because interchannel coupling is neglected outside r matc h. 

For the moderately anisotropic Mg + NH system stud- 
ied here, optimized MQDT with an interpolated Y ma- 
trix can provide numerical results in quantitative agree- 
ment with fully converged coupled-channel calculations. 
In future work, we will investigate the extension of this 
approach to more strongly coupled systems, with larger 
anisotropy of the interaction potential and more closed 
channels that produce scattering resonances. 
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